Introduction

The very nature of the global landscape in the 21st century can be characterized as one of uncertainty and unpredictability. Through the rise of globalization, a process compelled by investments and international trade, facilitated through the increase use of information technology, societies are becoming better integrated and more interconnected. However, one thing is certain: as globalization continues, strife amongst societies will persist. Our adversaries will continue to seek environment dominance and create an asymmetrical advantage by establishing regional hegemony to disrupt U.S. sanctuary at home or abroad. Bahrain is no exception as actors, both state and non-state, have sought to establish strategic control and influence over their population since 2011. These acts have resulted in violent and nonviolent protests occurring on a near daily basis within the kingdom. Therefore, developing a model to predict populace behavior is vital to U.S. strategic objectives.

Scenario

A team of analysts, LT Angeles, Capt Biles, and CPT Wren are assigned to 5th Fleet/US Navy Central Command (NAVCENT) located at NSA Bahrain. 5th Fleet has asked its analyst team to model the civil unrest within the island using open source data.

Collecting raw civil data within the Political, Military, Economic, Social, Information, Infrastructure, Physical Environment, and Time (PMESII-PT) variables allows the 5th Fleet analyst team to follow a commonly understood framework. By applying data mining techniques across open source websites, daily data beginning on January 1st, 2016, through December 31st, 2019, from across 191 PMESII-PT independent variables was gathered and collated. This time period resulted in a sample size, n, of 1461. By collecting data across PMESII-PT variables, one is able to develop a baseline understanding of an operational environment that is comprehended universally by military leaders. It is important to note that each operational environment is unique and that while some PMESII-PT variables are applicable, each environment requires its own data collection plan. While individual variables do not dominate every operational environment, they are all present and require careful consideration; ignoring one or more of these variables can negatively impact military missions and degrade a commander’s understanding of the operational environment.

Codebook

The predictors in master_df.csv utilized Political-Military-Economic-Social-Information-Infrastructure-Physical_Environment-Time (PMESII-PT) as a framework for data collection.

1. Political

• Political entities such as President Trump, Iranian President Hasan Rouhani, Bahrain’s Ministry of Interior, and Israeli Prime Minister Benjamin Netanyahu utilize various social media platforms to include Twitter. These data are captured using a Python script and a Twitter developer API and reported a daily count of the number of posts, favorites, and retweets.

2. Military

• Although arrest records are openly available via their department website, it is difficult to get historical data. As such, there is no data available in any dataframe.

3. Economic

• Data mining of Bahrain’s currency exchange rate reveals minimal fluctuation between 2.645 to 2.660 Bahraini Dinar to the United States Dollar. As a result, the impact of Bahrain’s exchange rate is minimal and excluded from this data set.
• Like most developed nations, Bahrain employs a stock market, the Bahrain Bourse or Bahrain Stock Exchange (BSE), where the issuing and trading of bonds, securities, and equities occurs through formal exchanges.
• In addition to the Bourse, the primary industries within Bahrain are: petroleum processing and refining, aluminum smelting, and iron palletization, which produce trackable indices and are tradable on the open market.
• Additional commodities such as gold, natural gas, and Bitcoin were also collected, which diversified the data collection portfolio.
• Some experts contend that one of the driving forces behind the Arab Spring in 2011 was the high cost of food. Therefore, the prices of corn, rice, wheat, cattle, soybeans, cocoa, cane sugar, and coffee were also tracked on the open market.
• When tracking these commodities on the open market, daily high, low, open, and close are all presented for trading days.
• However, on days when no trading occurred such as weekends and federal holidays, tracked prices were set to the previous day’s values. For example, Saturday and Sunday’s data are equal to Friday’s data.

4. No daily social or Infrastructure data collected as they are difficult to measure categorically and even harder to quantify.

5. Information

• According to the Arab Social Media Report 2015, mobile messaging application Whatsapp has a penetration rate of 85%, Facebook has a penetration rate of 81%, Instagram has a penetration rate of 51%, and Twitter has a penetration rate of 47% for Bahrain.
• Due to encryption protocols within Whatsapp, mining data and its usage are not possible.
• Various insurgent groups utilize social media as a medium for propaganda purposes to depict their activity and potentially serve as a recruiting tool.
• Some Twitter accounts post their messages in both English and Arabic while some only utilize Arabic.

6. Physical Environment

• The National Oceanic and Atmospheric Administration provides historical data for Bahrain at its international airport.
  
• Daily data collected included temperature, dew point, wind speed, gusts, and total rainfall.

7. Time

• Various aspects of time and the perception of them can affect and influence military operations both friendly and enemy within the environment.
  
• As each day passes, the numerical day, day of the week, and month are collected in passing.
  
• Holiday information is represented for both Bahraini political and religious holidays while moon illumination is produced by the U.S. Naval Observatory in Washington D.C. for each time zone around the world.

Our various target variables, 2016-01-01-2019-12-31-Bahrain.csv, are provided by The Armed Conflict Location & Event Data (ACLED) Project, which “is a disaggregated data collection, analysis, and crisis mapping project. ACLED collects the dates, actors, locations, fatalities, and modalities of all reported political violence and protest events across Africa, South Asia, Southeast Asia, the Middle East, Central Asia, and the Caucasus, Latin America and the Caribbean, and Southeastern and Eastern Europe and the Balkans. The ACLED team conducts analysis to describe, explore, and test conflict scenarios, and makes both data and analysis open for free use by the public.” Source: https://acleddata.com/about-acled/

For the initial stages of this analysis, the team will focus its efforts to predict the total number of riots, which is an aggregated sum of mob violence and violent demonstrations. Further analysis below conducts dimensionality reduction and attempting to predict that new target variable.

Exploratory Data Analysis (EDA): Part 1

The ACLED data set provided several different event and sub-event categories for study and analysis. (Notionally) 5th Fleet is interested in any event that could potentially cause a disruption in base operations. Analyzing the data can enable the team to identify patterns, discover anomalies, and generate hypotheses. Analyzing the raw data led the team to categorize social media posts, month, and weekday as discrete variables and holiday as binary. The team elects to conduct a more thorough EDA once it selects an event to model.

# Bar chart presenting event types over time along with their count
plot1 <- ggplot(by_date_type, aes(fill=event_type, y=count, x=event_date)) + 
    geom_bar(position="stack", stat="identity") + scale_x_date(date_breaks='year', date_labels='%y') + 
  labs(title='Bahrain Events (01-JAN-16 to 31-DEC-19)', x='Date', y='Count') +
  theme(legend.title=element_blank(), legend.position='right', legend.text=element_text(size=9))
  
ggplotly(plot1, tooltip=c('y', 'x', 'event_type')) 

From Plot1 we can see that protests make up the majority of our events. However, these protests can be divided into peaceful or violent and this warrants further analysis. Although violent and nonviolent protests occur on a near daily basis in the Kingdom of Bahrain, these protests are localized to predominately Shi’ite neighborhoods such as Diraz, Sitra, and Ma’ameer. Despite being predominantly Shi’ite districts, these neighborhoods present distinct social, economic, and political challenges when identifying and understanding the civil domain.

# Bar chart presenting sub event types over time along with their count
plot2 <- ggplot(by_date_type, aes(fill=sub_event_type, y=count, x=event_date)) +
  geom_bar(position="stack", stat="identity") + scale_x_date(date_breaks='years', date_labels='%y') +
  theme(legend.title=element_blank(), legend.position='bottom', legend.text=element_text(size=9)) + 
  labs(title='Bahrain Sub-events (01-JAN-16 to 31-DEC-19)', x='Date', y='Count')
  
ggplotly(plot2, tooltip=c('y', 'x', 'sub_event_type')) 

From Plot2 we can see that the majority of events are peaceful demonstrations. However, this conclusion leads to further discussion about peaceful protests turning violent or action caused by the government against civilians. The code analyzing this phenomenon is further below (EDA Question #3).

# GG plot for riots during our time frame
plot4 <- master_df %>%
  ggplot( aes(x=event_date, y=riot_count)) +
    geom_area(fill="#69b3a2", alpha=0.5) +
    geom_line(color="#69b3a2") +
    labs(title='Bahrain Riots (01-JAN-16 to 31-DEC-19)', x='Date', y='Count') + 
  scale_x_date(date_breaks='years', date_labels='%Y')
    
ggplotly(plot4, tooltip=c('y', 'x', 'event_type'))

In Plot4 we can see that although violent demonstrations occur nearly every day, there appears to be a cyclical nature for high volume events. Terrorists and insurgents often use a cycle to carry out their attacks. According to the Director of National Intelligence, the terrorism cycle roughly consists of re-attack surveillance, training, and rehearsals. Conversely, the spike in protests may be in response to government actions against civilians or simply pure coincidence. In the R chunk below we will explore the seasonality nature of these demonstrations (EDA Question #1). A time series is said to be stationary if it holds the following conditions true:

1. The mean value of time-series is constant over time, which implies, the trend component is nullified.
2. The variance does not increase over time.
3. Seasonality effect is minimal. Source: http://r-statistics.co/Time-Series-Analysis-With-R.html

In order to determine seasonality, we will utilize the Webel and Ollech (2019) test, which is an overall test for seasonality of a given time series. Source: https://cran.r-project.org/web/packages/seastests/seastests.pdf

# EDA Question #1: Is there evidence of a terrorist/insurgent attack cycle via seasonality analysis?

# Code sourced from https://cran.r-project.org/web/packages/seastests/vignettes/seastests-vignette.html

# Creates a dataframe for conversion into a time series object
tsData <- cbind(ts(master_df[, 1]), master_df[, 223])

# Creates a time series object
ts_obj <- ts_ts(ts_long(tsData))

# Performs a simple Webel and Ollech seasonality test
isSeasonal(ts_obj, test = 'wo', freq = 7)

# Checks for seasonality within a month
# Create a sequence of days from 1 to 30
day_seq <- seq(1,30, 1)
i <- 1 # Initialize a counter

# Create empty data containers for storage
kw_df <- data.frame() 

# Iterate over the length of 30 days
while(i <= length(day_seq)){
  # store the results of the Kruskal Wallace test
  kw_df[i,1] <- kw(ts_obj, freq = i, diff = T, residuals = F, autoarima = T)$Pval
  kw_df[i,2] <- i
  
  i = i + 1 # Increment the counter
}

# Change column names
names(kw_df)[1] <- "P-value"
names(kw_df)[2] <- "Day"

# Produce the trend and remainder data.
forecast::mstl(master_df$riot_count)

## [1] "Day  1  has p-value 1"
## [1] "Day  2  has p-value 0.568633129024278"
## [1] "Day  3  has p-value 0.0895051025727108"
## [1] "Day  4  has p-value 0.874688632633448"
## [1] "Day  5  has p-value 0.031469940685127"
## [1] "Day  6  has p-value 0.167023568005817"
## [1] "Day  7  has p-value 0.391981114540598"
## [1] "Day  8  has p-value 0.621282101342723"
## [1] "Day  9  has p-value 0.510489110739845"
## [1] "Day  10  has p-value 0.0849644787543207"
## [1] "Day  11  has p-value 0.159815642799602"
## [1] "Day  12  has p-value 0.516291197778634"
## [1] "Day  13  has p-value 0.0445863214864214"
## [1] "Day  14  has p-value 0.56186090500358"
## [1] "Day  15  has p-value 0.196047602054947"
## [1] "Day  16  has p-value 0.0558027515517882"
## [1] "Day  17  has p-value 0.953299557161738"
## [1] "Day  18  has p-value 0.657759469469912"
## [1] "Day  19  has p-value 0.378992512650113"
## [1] "Day  20  has p-value 0.212477444683822"
## [1] "Day  21  has p-value 0.818646678140778"
## [1] "Day  22  has p-value 0.23977688981302"
## [1] "Day  23  has p-value 0.570585359465244"
## [1] "Day  24  has p-value 0.308238021403119"
## [1] "Day  25  has p-value 0.123195116716367"
## [1] "Day  26  has p-value 0.0594705662306877"
## [1] "Day  27  has p-value 0.742004623410293"
## [1] "Day  28  has p-value 0.16403076361633"
## [1] "Day  29  has p-value 0.878928910906065"
## [1] "Day  30  has p-value 0.379182639938412"

## [1] "The estimated cycle is 5 days"

By utilizing the Kruskal-Wallace Test with 30 groups, we are able to determine if there are differences between two or more groups of an independent variable on a continuous or ordinal dependent variable (source: https://statistics.laerd.com/spss-tutorials/kruskal-wallis-h-test-using-spss-statistics.php). From the output above, we are able to conclude that there is sufficient evidence that riots occur on a 5 or 13-day cycle due to their low p-value. Further analysis of Plot5 contributes to this conclusion as it’s trending mean is approximately zero with an almost constant variance. As these attacks appear to be seasonal, the team lacks sufficient evidence to pinpoint the exact duration of the insurgent planning cycle.

Speculation amongst the analyst team at 5th Fleet (notional) believes that these protests are localized within predominately Shi’a neighborhoods who believe they are oppressed by the government. Identifying the towns with the most disruptions will enable effective force protection measures and ensuring minimal risk to the base.

EDA Question #2: Where are the protests occurring and what is the subsequent risk to the base and personnel?

# EDA Question #2: Where are the protests occurring and what is the risk to the base?

# Aggregate protests by location and exclude non-violent protests
by_location_type <- raw_df %>% group_by(location, sub_event_type != 'Peaceful protest') %>% dplyr::summarise(count = n())

# Change column name to Boolean_Violence
names(by_location_type)[2] <- "Boolean_Violence"

# Filter data that only accounts for violent protests only
filtered_location_df <- data.frame(by_location_type[by_location_type$Boolean_Violence == TRUE,])

# Sort the dataframe above with the highest number of events
filtered_location_df <- head(filtered_location_df[order(filtered_location_df$count, decreasing = TRUE),],10)

# Top 10 locations for violent/nonviolent protests
filtered_location_df['longitude'] <- NA
filtered_location_df['latitude'] <- NA
violent_locs <- filtered_location_df %>% pull(location)
for (loc in violent_locs) {
  filtered_location_df[filtered_location_df$location==loc,]['longitude'] <- raw_df[raw_df$location==loc,]$longitude[1]
  filtered_location_df[filtered_location_df$location==loc,]['latitude'] <- raw_df[raw_df$location==loc,]$latitude[1]
}

# All locations for violent/nonviolent protests
by_loc_type_df <- as.data.frame(by_location_type)
unique_locs <- by_loc_type_df %>% pull(location)
unique_locs <- unique(unique_locs)
by_loc_type_df['longitude'] <- NA
by_loc_type_df['latitude'] <- NA
for (loc in unique_locs) {
  by_loc_type_df[by_loc_type_df$location==loc,]['longitude'] <- raw_df[raw_df$location==loc,]$longitude[1]
  by_loc_type_df[by_loc_type_df$location==loc,]['latitude'] <- raw_df[raw_df$location==loc,]$latitude[1]
}

# Bar plot for top 10 violent demonstrations
plot6 <- ggplot(filtered_location_df, aes(x=reorder(location, -count), y=count, fill=location)) + geom_bar(stat='identity') + labs(title='Top 10 Violent Demonstrations (01-JAN-16 to 31-DEC-19)', x='Location', y='Count') + geom_text(aes(label=count), vjust=1.6, color="white", size=3) + scale_fill_brewer(palette='Paired') + theme(legend.position='none', axis.text.x=element_text(angle=45, hjust = 1))

plot6

Sys.setenv('MAPBOX_TOKEN' = 'pk.eyJ1IjoiZGF2aWR3cmVuIiwiYSI6ImNrYWxqNG8xZTB0bW0ycmxiNGI2NW9jenAifQ.XspCPRtcou2DCia28IF3og')

# Bubble map for where violent demonstrations are occurring
plot7 <- plot_mapbox(maps::canada.cities) %>%
  add_markers(data= by_loc_type_df[by_loc_type_df$Boolean_Violence==TRUE,],
    x = ~longitude,
    y = ~latitude, split=~location,
    size=~count,
    color = ~count,
    text=~paste('Count:', count)) %>%
  layout(
    showlegend=FALSE, title='Bahrain Riots (01-JAN-16 to 31-DEC-19)',
    mapbox=list(center=list(lat=median(raw_df$latitude), lon=median(raw_df$longitude)), zoom=9))

plot7

# Bubble map for where non-violent demonstrations are occurring
plot8 <- plot_mapbox(maps::canada.cities) %>%
  add_trace(data= by_loc_type_df[by_loc_type_df$Boolean_Violence==FALSE,], type='scattergeo',
    x = ~longitude,
    y = ~latitude, split=~location,
    size=~count,
    color = ~count,
    text=~paste('Count:', count)) %>%
  layout(
    showlegend=FALSE, title='Non-violent Demonstrations (01-JAN-17 to 31-DEC-19)',
    mapbox=list(center=list(lat=median(raw_df$latitude), lon=median(raw_df$longitude)), zoom=9))

plot8

The bar graph compares the number of riots of the top ten cities in Bahrain that have the highest count; Nuwaidrat experienced the most violent demonstrations at 387. From the bubble map of violent demonstrations, we see that not very many protests actually occur around NSA Bahrain; there were only 91 violent protests northwest of NSA Bahrain in Juffair. Most of the protests occurred in the south-southwest region of Bahrain, largely focused on the cities of Sitra and Ma’ameer. Conversely, large numbers of peaceful protests are more spread out through the country. While Ma’ameer still experienced many peaceful protests, cities like Budaiya, Jid Hafs, and A’ali have also experienced a large amount. Additionally, we see that not very many peaceful protests occurred around NSA Bahrain, with the closest city being Juffair at 25 protests.

# EDA Question #3: What is the relationship between nonviolent protests and violent demonstrations?

# This code determines the correlation between peaceful protests and violent demonstrations
cor.test(master_df$riot_count, master_df$Peaceful_protest_count)

## 
##  Pearson's product-moment correlation
## 
## data:  master_df$riot_count and master_df$Peaceful_protest_count
## t = 21.132, df = 1459, p-value < 2.2e-16
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  0.4438215 0.5224032
## sample estimates:
##       cor 
## 0.4840878

# Create a dataframe for correlation plotting
corr_df <- data.frame(master_df$riot_count, master_df$Peaceful_protest_count)

# Change column names
names(corr_df)[1] <- "Riots"
names(corr_df)[2] <- "Peaceful_Protests"

# Code sourced from: http://www.sthda.com/english/wiki/correlation-test-between-two-variables-in-r
# Create a plot
plot9 <- ggscatter(corr_df, x = "Riots", y = "Peaceful_Protests", 
          add = "reg.line", conf.int = TRUE, 
          cor.coef = TRUE, cor.method = "pearson",
          xlab = "Violent Demonstration Count", ylab = "Peaceful Protest Count", title = 'Violent and Peaceful Demonstrations')

plot9

## `geom_smooth()` using formula 'y ~ x'

# Shapiro-Wilk normality test for peaceful protests
shapiro.test(corr_df$Peaceful_Protests)

## 
##  Shapiro-Wilk normality test
## 
## data:  corr_df$Peaceful_Protests
## W = 0.72153, p-value < 2.2e-16

# Shapiro-Wilk normality test for violent demonstrations
shapiro.test(corr_df$Riots)

## 
##  Shapiro-Wilk normality test
## 
## data:  corr_df$Riots
## W = 0.71688, p-value < 2.2e-16

We see the correlation calculated above is 0.4840878, which is not above the widely accepted standard of 0.7 or -0.7. However, given the nature of the events measured, we see that the correlation coefficient is insightful: not every non-violent protest will become violent. Human actions from both parties and other external influences can either escalate or de-escalate a non-violent protest.

The team scrutinized a subset of the data and looked at individual events and at this point believes that non-violent events and violent demonstrations are not double-counted implying they are IID. Furthermore, the team believes that both peaceful protests and violent demonstrations are both normal via the Shapiro-Wilk normality test. As a result, the team believes that incorporating non-violent events may serve as a useful predictor variable.

From the data above we can see that battles and attacks do occur between Ministry of Interior (MOI) forces and the civilian population. The late Martin Luther King Jr stated that “violence begets violence.” Research into this topic leads the team to its fourth EDA question: Does violence beget violence?

# EDA Question 4: Does violence beget violence?

# Compute the correlation between violence against civilians and violent demonstrations
cor.test(master_df$riot_count, master_df$VAC_count)

## 
##  Pearson's product-moment correlation
## 
## data:  master_df$riot_count and master_df$VAC_count
## t = 0.97831, df = 1459, p-value = 0.3281
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  -0.02571446  0.07678795
## sample estimates:
##        cor 
## 0.02560404

From the correlation coefficient, we can see there is no discernible relationship between violence against civilians and violent demonstrations on the same day. However, additional analysis will need to occur to determine whether or not violence against civilians serves has an impact in the coming days as news of an attack would spread amongst the population and incite emotional responses.